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We consider a fluid of hard spheres bearing one or two uniform circular adhesive 
patches, distributed so as not to overlap. Two spheres interact via a "sticky" Baxter 
potential if the line joining the centers of the two spheres intersects a patch on each 
sphere, and via a hard sphere potential otherwise. We analyze the location of the 
fluid-fluid transition and of the percolation line as a function of the size of the patch 
(the fractional coverage of the sphere's surface) and of the number of patches within 
a virial expansion up to third order and within the first two terms (CO and CI) 
of a class of closures Cn hinging on a density expansion of the direct correlation 
function. We find that the locations of the two lines depend sensitively on both the 
total adhesive coverage and its distribution. The treatment is almost fully analytical 
within the chosen approximate theory. We test our findings by means of specialized 
Monte Carlo (MC) simulations and find the main qualitative features of the critical 
behaviour to be well captured in spite of the low density perturbative nature of 
the closure. The introduction of anisotropic attractions into a model suspension 
of spherical particles is a first step towards a more realistic description of globular 
proteins in solution. 
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I. INTRODUCTION 

The idea of modeling fluids as systems of spherical particles with orient at ionally depen- 
dent attraction dates back at least as far as Boltzmann, who envisaged chemical attraction 
between atoms only when "their sensitive regions are in contact".- Models of this type, 
featuring patchy interactions, are currently experiencing renewed relevance in the context 
of colloidal and biological matter in contrast to their original conception in connection with 
fluids of atoms and small molecules.-^'^^'^'^-^^'^^'^^^^^'^^^ 

The new interest arises for various reasons. On the technological side, patchy particles 
give the possibility of designing self-assembling nano-scale devices through anisotropic dec- 
orations of the particle surface by means of organic or biological molecules.— 'i^'i^ Nature 
provides inspiration for what might be achieved in this area, a particularly elegant exam- 
ple being the self assembly of virus capsids. These protein shells are monodisperse and 
highly symmetric, and are composed of identical subunits. Simplified descriptions of icosa- 
hedral virus capsids are currently being formulated using spherical subunits with directional 
interactions,— and the possibility of adopting similar schemes to self assemble other tar- 
get structures is being explored.— This level of organization inevitably requires a certain 
specificity in the interactions between the subunits as well as measures to prevent further 
aggregation of the assembled objects. 

Less specific patchy interactions give rise to associating fluids containing a distribution 
of cluster sizes or extended gel-like networks. The key feature of such systems is a set 
of point-like sites on the particle surface, leading to strongly directional bonding with a 
maximum of one bond per site.-'^ii^'iiii^'i^ This type of interaction has proven invaluable 
in elucidating the interplay between fluid-fluid and sol-gel transitions. One advantage of 
these models is that powerful analytical tools are available for them, such as Wertheim's 
thermodynamic perturbation theory,— which yields accurate results under experimentally 
realistic conditions.—"^ 

In contrast to these models with attractive spots, one can envisage particles that interact 
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through larger attractive regions on their surface, for example, globular proteins with patches 
of hydrophobic (nonpolar) amino acids exposed at the surface. Isotropic potentials have 
been remarkably successful in modeling the phase diagrams of certain proteins,— "^"^ but 
it seems that not all features of their coexistence curves can be properly explained by such 
simple interactions.— In this sort of system, it seems more appropriate to consider regions 
with short-range attractive forces^'^'^'^'^'^ii^ rather than site-site bonds. These attractive 
patches are capable of sustaining as many "bonds" as permitted by geometry. The size of 
the patch therefore becomes an important new parameter that does not arise in most work 
on associating fluids. We note, however, that for sufficiently narrow patches, the two models 
become essentially equivalent. 

In the present work, we focus on a simple yet physical model which is a variation of those 
treated in.— We consider uniform circular patches distributed on the surface of the sphere 
in such a way that they do not overlap. The patches are delimited by circles which can 
be defined by the associated solid angles. Two particles experience an adhesive attraction 
only when a patch on one sphere touches a patch on the other. The adhesion is of Baxter's 
type,— i.e., the attraction has infinitesimal range, acting only when the particles are exactly 
in contact, as described in the next section. This model has the advantage that it can be 
tackled with analytical tools, unlike most other models for which not even the isotropic 
analogue bears this appealing feature. Various issues arise in this sort of model relating to 
the stability of the liquid phase with respect to crystalline solid phases, and these points 



have been studied in Refs. 
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The integral equation theory of fiuids with an angularly dependent pair potential is 
complicated by the fact that the pair distribution function is also angularly dependent.— In 
the general case one must appeal to the symmetries of the fiuid (translational invariance, 
rotational invariance, invariance under permutation of like particles, and invariance under 
the symmetry operations of the individual particles and of the correlation functions) in 
order to simplify the problem.— "^"^ In some cases, it is possible to factorize the angular 



dependence of the Ornstein-Zernike (OZ) equation. For examp 
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e, the factorization for a 
it is shown how to solve 



fiuid of dipolar particles has long been known,— and in Ref. 
the dipole-dipole angular distribution of attraction in the adhesive limit within the Percus- 
Yevick (PY) framework. However, the dipolar case hinges on exploiting a special property 
of this particular angular distribution that is particularly useful for the angular convolution 
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in the OZ equation. In contrast, for an angular dependence with discontinuities, such as 
the circular patches treated here, any approach relying on a spherical harmonic expansion 
would prove a formidable task due to the large number of terms necessary to capture the 
discontinuities. 

In the present model we therefore follow a different route based on two parallel and re- 
lated schemes. We first perform a virial expansion up to the third virial coefficient. We then 
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proceed to study a class of closures (denoted CO, CI,...) which were proposed in Ref. 
and are based on a density expansion of the direct correlation function. In particular, the 
zeroth-order term (CO) turns out to be equivalent to a modified mean spherical approxima- 
tion, whereas the first-order (CI) is known to provide the correct third virial coefficient.— 
Within both schemes we study the thermodynamics, radial distribution function and perco- 
lation threshold, and compare with specialized Monte Carlo simulations which were recently 
devised to this aim.— By varying the size of the adhesive patches and by selecting between 
one patch and two diametrically opposite patches, we are able to investigate the roles of 
both the total surface coverage and the geometrical distribution of the adhesion. In both 
the one- and two-patch cases we can change smoothly between small sticky spots, capable 
of making only one bond each, and the isotropic adhesive sphere. We find that the position 
of the critical fluid-fluid transition line and the percolation threshold are both sensitive to 
the surface coverage. At fixed coverage, there is also a dependence on the way in which this 
adhesion is distributed. 

Our results can be compared and contrasted with the recent work of Bianchi et al.,— 
who consider the maximum number of bonds per particle, rather than the fractional surface 
coverage, as the key parameter controlling the location of the critical point. In the present 
work we are able to tune both effects, thus illuminating their specific roles in the location 
of critical points. 

The remainder of the paper is organized as follows. In section |Tl] we introduce the model 
while section Hm contains a description of the analytical and numerical used tools. Results for 
the radial distribution function, fiuid-fiuid transition and percolation threshold are included 
in sections llVt IVl and IVIt respectively. Finally, in section IVIII the inclusion of an adhesive 
background is discussed, and conclusions and an outlook are contained in section IVIIII 
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II. DEFINITION OF THE MODEL 



A. Baxter model with orientationally dependent adhesion 



We start with some general remarks on the orientational dependence of a three- 
dimensional homogeneous fluid of hard spheres with adhesive pairwise interactions. Let 
Ti be the coordinates of the i^^ particle {i = 1,2,3,.. .) and assume that the patch distribu- 
tion on the sphere has cylindrical symmetry so that its orientation in space is determined by 
a unit vector Sj rigidly attached to it. Then Sj = (sin 6i cos ipi, sin 6i sin ipi, cos 6i) where 6i and 
(fi are the polar and the azimuthal angles with respect to a fixed reference frame (see Fig. [1]). 
As usual, we introduce the relative coordinates ri2 = r2 — ri and the associated distance 
^12 = |ri2|; and work with the following short-hand notation: (1, 2) = (ri2, ^i, (fi, 62, ^^2) and 
Qi = {6i,ipi) for the orientation of Sj. The orientation of = rij/vij with respect to the 
same frame of reference will be denoted by Qij = {6ij, (pij). 

The particles interact through a pair potential 0(1, 2), defined below, which is a general- 
ization of Baxter's sticky hard sphere (SHS) limit^ to orientationally dependent interactions. 
We start with 
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where f3 = l/lksT) {ks being Boltzmann's constant and T the temperature), a is the 
diameter of the spheres, and e(l,2)/r is a dimensionless adhesion coefficient. We define 
0(1, 2) through the following limit on the Boltzmann factor e, 

e(l,2) a 



e(l, 2) = exp[-/30(l, 2)] = lim exp[-/3$(l, 2)] = e(ri2 - a) 



-5(ri2 - a) . (2) 



where G(-) is the Heaviside step function and 5{-) the Dirac delta function. When e(l, 2) = 1 
we recover the usual Baxter SHS model and hence the only orientational dependence is 
included in the definition of e(l, 2). It is easy to see that e(l, 2) cannot be a simple function 
of §1 and §2 but must also include a dependence on V12 = ri2/ri2, in order to avoid a trivial 
corresponding states rescaling. This point is discussed in Appendix El In the present work 
we shall address a type of orientational dependence which was introduced by Kern and 
Frenkel^ following a previous suggestion by Jackson et al- 
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B. Patchy sticky hard spheres 

Consider a single hard sphere having one or more identical adhesive circular patches 
distributed on its surface in such a way that they do not overlap with one another. The 
size of the patch can be specified by the angular amplitude 25 as shown in Fig. O The unit 
vector sf"* identifies the direction from the center of particle i to the center of patch p on the 
surface (p = 1, . . . , n, the total number of patches). The sticky area is then given by points 
f on the surface of the particle such that the angle between s^-^^ and f is smaller than 6. 

In conjunction with Eq. ([2]), the adhesive part of the interaction between two particles 
acts only if their point of contact lies inside a patch on each particle, as depicted in Fig. [3] 
for the case of a single patch (n = 1). Therefore, e(l, 2) = e(si, S2, ri2) can be written as 

{1 if s^^^^ ■ ri2 > cos 5 and — s^^'^^ ■ Yio > cos 5 for some combination (pi, no) 
otherwise. 

Each patch occupies a portion of the sphere's surface covered by the solid angle 27r(l — 
cos 6) and a fundamental role will be played in our discussion by the fraction of solid angle 
(i.e., the coverage) associated with 6, namely 

Xo(5) = ^(l-cos5) = sin2('^V (4) 



2 V2. 

III. ANALYSIS OF THE MODEL 

A. Analytical solution 

We now tackle the analytical solution of this problem based on two simple approximations: 
the virial expansion and the Cn class of closures. 

1. Virial expansion 
As shown in Appendix B, the first two virial coefficients for this model are 

b2 = B^/vo = 4 - 12-^ , (5) 
63 = B,/vl = 10 - 60|L + 144^ - 96 . (6) 
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where Vq = ira^/G is the volume of a sphere and 




(7) 
(8) 
(9) 



where we have defined the angular average (with dfl = dfl/An) 
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(10) 



Here, 6ij is the angle between fjj and (which can be chosen along the z axis), and 
is always associated with a delta function that forces spheres i and j to be in contact. Note 
that in Eq. ([6]) the effect of anisotropy is embedded in the parameters Xi)X2,X3 defined in 
Eqs. (171), (JHl), and ([9]), and that these parameters are therefore functions of 6 and n. The 
isotropic case is recovered when all x's equal 1. We remark that the expression for X2 involves 
an average over the relative orientations while there is an overlap between spheres 1 and 
2, each of which is simultaneously in contact with sphere 3. Under such conditions there is 
always a maximum possible angle vr/S for ^13 and this gives rise to the normalization factor 
of 4 in Eq. (El). 

If one limits the expansion to the second virial coefficient, a law of corresponding states 
based on the rescaling r t/xi between the patchy and the isotropic SHS models holds 
true. This correspondence breaks down even at the level of the third virial. 

It is easy to see that xi = ""-^Xo ^^is is simply the product of the separate coverages 
on each sphere. A calculation of X2 and xs is much more laborious and can be found in 
Appendix O for the case of a single patch. The final result in this case is 



where the coefficients Qi and Ri are given in Appendix O For 6 > bn/G it is possible 
to have three mutually bonded spheres with the patch vectors pointing either inward and 
outward (see Fig. H]). Note that for the isotropic limit 5 = tt all (^ = 1? 2, 3) are equal to 
1 as they should be. The three Xi coefficients are plotted in Fig. [5l as functions of 6. 



Xi = xl , 
X2 = xlQii^) , 
X3 = RliS) , 



(11) 

(12) 
(13) 
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For spheres with two diametrically opposite patches, each of width 5, one finds 

Xi = 4xo , (14) 
X2 = ^xlQ2i5) , (15) 

X3 = RUS) , (16) 

where the coefficients Q2 and R2 are given in Appendix [Dl Note that in this case when 
(5 > 7r/3 it is also possible to have sphere 1 in contact with spheres 2 and 3 through different 
patches as shown in Fig. [6l The three Xi coefficients are plotted in Fig. [7] as functions of 6. 
The virial expansion of the excess free energy density is 

Prvo = hv' + Ihr]-' + . . . , (17) 

where rj = pvo is the hard sphere packing fraction. This allows the calculation of the 
corresponding pressure and chemical potential, 

f3P{T, ri)vo = 7] + h27f + hrf + . . . , 

3 

/3/i(r, 7]) = ln(AVt;o) + \nr] + 2b2V + i^hv^ + • ■ ■ , 
where A is the de Broglie wavelength. 



2. Integral equations within the Cn closures 

While the virial expansion only allows a limited low-density region of the phase diagram 
to be probed, the integral equation approach is much more powerful in this respect. The 
trade-off is, of course, that since the OZ equation involves the total correlation function h and 
direct correlation function c, both of which are unknown, it can be solved only after adding 
a closure, that is a second, approximate, relationship involving /i, c, and the pair potential. 
In this section we discuss a particular class of these closures (denoted Cn hereafter) which 
have already been exploited in the isotropic case and have proven to provide reasonably 
good predictions even for intermediate densities.— 

The OZ equation for a homogeneous ffuid of molecules interacting through anisotropic 
pair potentials is 

/i(l, 2) = c(l, 2) + p / rf(3)c(l, 3)/i(3, 2) , (18) 
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where d{i) = dvidfli. More explicitly [see Eq. (fTOj) ] 

2) = c(l, 2) + p y rfr3(c(l, 3)/i(3, 2))f,3 . (19) 

In a homogeneous fluid, translational invariance of any correlation function implies that 
one can introduce reduced coordinates = r2 — ri and ri3 = ra — ri 

/i(ri2,fii,fi2) = c(ri2,fii,fi2) +P J rfr3(c(ri3,fii,fi3)/i(r32,ri3,^2))n3 • (20) 

The presence of the convolution makes it convenient to Fourier transform this equation with 
respect to the position variable r. This yields for the corresponding functions (indicated 
with a hat) in Fourier space k 

h{k, Q,, Q2) = a(k, fil, ^2) + p(c(k, fil, Q,)h{k, Q,, Q2))ns ■ (21) 

The additional complication with respect to the isotropic case is the presence of the orien- 
tational average of the product appearing in Eq. fl2Tl) . In order to make progress, we use a 
simple angular decoupling approximation 

{c{k,Q,,Q;)h{k,Q,,Q2))n, ^ {c{k,Q,,Q3))nXh{K^3,^2))ns ■ (22) 



As discussed in Ref. 
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Cn closures are based on a density expansion of the cavity function 
y{l, 2) which is related to the radial distribution function g{l, 2) by g{l, 2) = y{l, 2)e(l, 2). 

In the presence of anisotropy, all correlation function clearly depend upon the solid angles 
Qi, Q2, ^12- It is then customary to consider— the corresponding angular averaged quantities 
g{r = ru) = {g{l, 2))^^^^^^^^^ and similarly for y{r = r^). Within Cn closures, for r > a, the 
radial distribution function g{r) coincides with the cavity function y{r). A density expansion 
of the cavity function yields 

y{r) = l+py,{r) + ... , (23) 

where 

l/i(ri2) = I drs (/(l,3)/(3,2))^ (24) 

Calculation of Eq. ( l24l) proceeds using arguments akin to those presented in Appendix [HI 
which are based on the decomposition in Eq. (1B2I) . The integral in Eq. (IMl) then splits into 
three integrals containing the various combinations of the HS and the sticky parts of the 
Mayer function as in Eq. flB2p . 
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B. Monte Carlo algorithms for Baxter-like potentials 

Monte Carlo simulations of adhesive hard spheres require particular care even in the 
isotropic case because of the singular nature of the potential. For completeness we summarize 



the main ideas below, deferring to Ref. |37| for the details. 

Conventional Monte Carlo displacements of a SHS would fail because the bonded states 
between particles occupy an infinitesimal volume of configuration space (and so would never 
be located by random displacements) but have infinite strength (and so would never be 
broken). The solution is to compare the integrated weights of the various bonded and un- 
bonded states, which are finite. Specialized algorithms that exploit this approach have been 
devised for the canonical ensemble^*^ and were subsequently extended to the grand canon- 
ical ensemble.—*^ The latter is particularly convenient for identifying the critical point.— 

The Monte Carlo algorithm for isotropic adhesive spheres can be modified to deal with the 
patchy case by incorporating the anisotropy in the acceptance criterion for trial moves. Trial 



moves are attempted as described in detail in Ref. |41| as though the spheres were uniformly 
adhesive. Once the trial position of the displaced particle has been chosen, a uniformly 
distributed random orientation is selected. The move is then accepted only if an overlap 
of hard cores is not generated (as in the isotropic case) and if all contacts specified in the 
trial configuration have patches suitably aligned to make the required bonds. This scheme 
produces the desired Boltzmann distributional and is applicable to an arbitrary arrangement 
of patches. However, it becomes inefficient when the total adhesive coverage of the sphere 
is small because the random generation of orientations is then unlikely to lead to patches 
being aligned with bonds, leading to a high rejection rate. 

IV. STRUCTURE 

In the following we shall compare predictions from the combined Cl-orientational mean 
field approximation and virial expansion with the results of Monte Carlo simulations. 
One finds that yi{r) is different from zero only in the region < r < 2a and 
^ /r \ /r \2 /r \ X2{^) ^ 



yJr) = — {- + A](--2] +^27r --2 +-^^27r-, (25) 

12 V a J \a J 12t Ka J (12r)2 r' ^ ^ 

where 

Mr) = (e(l,2)e(l,3))^^f^^f^j0,,=2arcsin(r/2a). (26) 
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In order to compute the angularly averaged radial distribution function g{r) we have 
solved the full OZ equation fl20|) within the CI closure supplemented with the decoupling 
approximation (l22l) . The Wertheim-Baxter method^ and the CI closure combine such that 
only the cavity function at contact depends upon the angular coefficients XiyX2,X3- The 
solution for Baxter function is 



where 



q(r) = [a{r' -a')/2 + ba{r-a) + q^a']e{a-r) r>0 (27) 



2(1 - r/)2 1 - r/ ' 

yCl 

^ci = ^^ci (^^^ ^ 2))n„n2,ni2 = + • (31) 

The coefficients yo,i = 2/o,i(t) are related to the reduced virial coefficients 62,3 by Eq. fj^Tj) . 
below. Therefore we can read off their values 

yo = Xi , (32) 

In order to extract the numerical behavior of the radial distribution function g{r) we have 
employed a discretization method due to Perram^ to compute the numerical integral 



rh{r) = —q' [r) + 2ti p / du q{u){r — u)h{\r — u\) . (34) 

Jo 

For the one-patch spheres the result is reported in Fig. [8]for various values of 5, at r/ = 0.4 
and r = 0.2. The choices 5 = and 5 = correspond to the limiting cases of pure HS and 
isotropic SHS, respectively. Upon decreasing the size of the patch, the behavior smoothly 
interpolates between these two expected. The characteristic jump in g{r) air = 2a 

in the isotropic SHS model^ can be explicitly computed within the CI integral equation 
closure to be 



g{2a+)-g{2a~) = -Qr^[f^/{l2T)f 



(35) 
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The jump is also present for intermediate values of S and gradually fades out toward the HS 
result as illustrated in Fig. [HI In order to assess the precision of the predictions of the Cn 
closures, in Fig. [9] we compare the radial distribution generated by both CO and CI closures 
with Monte Carlo simulations (right-hand panel) and with the corresponding isotropic case 
(left-hand panel). The CI closure is seen to follow the Monte Carlo behavior well over the 
range of the ratio r/a considered, in both the isotropic and anisotropic cases. 

Like the Percus-Yevick solution of the isotropic model, the CO and CI closures fail to 
capture certain 5-function and step discontinuities in the radial distribution function,— such 
as those visible in the Monte Carlo results in the range a < r < 2a in Fig. [9l These features 
arise from clusters in which the distance between two particles is fixed or limited indirectly 
by a sequence of adhesive bonds, for example, the outermost pair of particles in face-sharing 
tetrahedra. These clusters are sampled correctly by the Monte Carlo simulations.— 

V. FLUID-FLUID COEXISTENCE CURVE 

An interesting issue, both from the theoretical view point and for the possible implica- 
tions in predictions of experimental phase transitions in solutions of globular proteins, is the 
determination of the fiuid-fiuid coexistence curve for the patchy sticky hard spheres, which 
we now address. We note that the dense fluid (liquid) phase, though often long-lived, is gen- 
erally only metastable for systems of particles interacting through sufficiently short-ranged 
isotropic attractive forces. However, it has recently been predicted that a thermodynami- 
cally stable liquid will be recovered if the coordination number of the particles is restricted to 
a maximum of six or less.— Nevertheless, it seems that the specific details of the interactions 
must be taken into account before a firm conclusion can be drawn for a particular model.— 

A. Virial expansion 

In order to find the coexistence or binodal line we need to solve for ?7i(r) and ?72(t) the 
following set of equations. 



P(r,r/i) = P(r,r/2) , 
/i(r,r/i) = fi{T,r]2) ■ 
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A straightforward use of Eq. (fTTl) for the excess free energy density to this aim, however, 
yields meaningless results even at moderate densities as one could have expected from the 



outset. A way out of this problem was proposed in Ref. |47| in the context of polydisperse 
SHS fluids. The idea hinges on a modification of the Carnahan-Starling expression for the 
HS excess free energy density,— 



= T^^l^r,^ , (36) 



so that it matches the patchy SHS result up to third order in density. A possible choice is 

2 2 

Prvo = (c - l)r/ln(l -rj) + + Cj^^^ , (37) 

where c and d are parameters to be determined by expanding to third order in density and 
matching to Eq. (fT7I) . We then find 

h -262 + 1 



3 

&2-I 



(38) 



d = . (39) 



The pressure and the chemical potential are then 
2I + 3d{l - r]) + ri[r] - 2 + c(3 - r])] 



PPVo = 7] + 7] 



wa3/ N 1 / N (1 + c + 6d)r/ - (2 - 2c + 9rf)r/2 + (1 - c + 3c?)r/3 
^/i = ln(A7t;o) + ln?7 + (c- l)ln(l -77) + ^ — ^ Ji \ /i 



respectively. In the limit c = d = 1 Eq. fl?r|) reduces to Eq. flS^ as expected. 

The behavior of the binodal line as a function of 6 is shown in Fig. [TOl As 6 decreases, the 
coexistence region shrinks as expected, since HS fluids {6 = 0) admit only a single phase.— 



B. CI integral equation 

An alternative route is to start form the excess free energy stemming from the energy 
route of the CI approximation^ 

prvo = Prjvo - (62 - 6f )r/^ + (63 - &f )y , (40) 

where b^^ = hnij 00). The rescaled virial coefficients 6.„ = 6„(r) can be related to the 
values of the corresponding coefficients of the expansion for the cavity function at contact 
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y = Vo + ViV + + • • • by means of the relation (see e.g., Ref. l36| and references therein): 

n — 1 (IT 

Hence, we have 

PPvo = ri + 2ri — r] + [ -5— + ^ - — ^ V , (42) 



(?7 — 1)^ r V T 18 r 



= ln(AV.o) + Inr, + ^^^3^ " y2r/ + ^-5^ + ^ ' ) ± ' (^3) 

We remark that these results need no further orientational approximations as all effects of 
anisotropy are exactly included in Xi,2,3- 

The results for the binodals are shown in Fig. [TOl For each of the adhesive coverages 
depicted, both theoretical treatments predict that the binodal line for the two-patch model 
lies above its counterpart for a single patch. This could be expected on physical grounds, 
since a more distributed region of adhesion usually facilitates the aggregation process. A 
closer analysis, however, indicates that this is not always the case. This is shown in Fig. [TT] 
where we report the change in the critical point as a function of the adhesive coverage nxo 
of the sphere surface. The difference between one and two patches decreases as 6 increases, 
as expected, but it is clearly visible through the whole range of existence. Remarkably, 
there is an inversion of the two curves around 80% coverage. For approximately nxo > 0.8 
the critical temperature for two patches lies above its one-patch counterpart. This means 
that the fluid-fluid transition line is encountered at higher temperature when an identical 
adhesive coverage is distributed over two spots rather than consolidated into a single big 
patch. However, this is no longer true when the size of the patch becomes too small. The 
reason is that under such conditions, it is then possible to bind three or more particles within 
a single patch, whereas at most two particles (one on each of the two opposite patches) can 
be attached in the two-patch case. 

An additional noteworthy feature of Fig. [11] is the existence of a kink at 5 = 7r/3 (50% 
coverage) in the two-patch case. The kink is related to the degeneracy illustrated in Fig. [61 
A similar kink occurs in the single-patch curve aX 6 = bir/G (93.3% coverage). Again, this 
stems from a degeneracy, as illustrated in Fig. [H 

Due to the inherent difficulty of tracing out critical temperatures to low coverage, we 
have been unable to verify the crossing of the one- and two-patch critical temperatures by 
Monte Carlo simulations. Fig. [11] reports MC results down to around 60% coverage and the 
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two-patch critical temperature is always above the one-patch case. We suspect, however, 
that an inversion might still occur in MC simulations, but at lower values of Xo? as yet 
inaccessible to our simulations. 

The sensitivity of the shape and location of the coexistence curve to the geometry of 
the adhesive distribution is quite a remarkable feature of this archetypal patchy model. It 
seems likely, therefore, that a proper understanding of experimentally determined phase 
diagrams of globular proteins should take into account the nonuniformity of their surfaces 
and consequently of their interactions. 



VI. PERCOLATION 



A further interesting issue, already discussed in the context of the isotropic model, is the 
percolation threshold,—"^ to which we now turn. 



A. Virial expansion 

In one-patch systems only dimers can form for 6 up to n/Q, while clusters of any size 
are in principle possible above this threshold. In order to study the percolation threshold 



we can use the strategy devised in Ref. 151|- Based on the definition of the connectedness 
correlation function (see below), the percolation threshold is signaled by the divergence of 
the mean cluster size 

S=l + pj dv,, (/i+(ri2, fii, ^2))^^^^^ , (44) 

where h'^ is the pair connectedness function, which is related to the direct connectedness 
function by the Ornstein-Zernike equation. Both are related to the connected part of the 
Mayer function /+(1,2) = /(I, 2) - /hs(1,2) as given in Eq. ({62]), /hs(1,2) being the HS 
part as given by Eq. ( 1B3I) . As in the case of Eq. (l22l) we assume that 

(c+(k,fii,fi3)/i+(k,fi3,^^2))n3 - (c+(k,fii,fi3))n3(/i+(k,f^3,^^2))n3 • (45) 

The average Fourier transform of the direct connectedness function c^{k) = 
{c^{k, Qi, fi2))ni,n2 a.t k = then identifies the threshold by the equation 

pc+(0) = 1 . (46) 
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Upon power expansion in the density we have 

oo 

c^O) = E (0)P"-^ = ct (0) + pet (0) + O (p2) . (47) 



n=2 



Using the above decomposition of the Mayer function the first two coefficients are found to 
be 

4(0) = /rfri,(r(l,2)>^^_^^, (48) 

£3+ (0) = j dv^^dvr,{ [/+(1, 2)/hs(1, 3)/hs(2, 3) + 3/hs(1, 2)/+(1, 3)/+(2, 3) + 

/+(l,2)/+(l,3)r(2,3)]>^^^^^^^^. (49) 
An analysis following that in Appendix [B] then yields: 

c+(0)/^o = 24^ , (50) 

ctm/vl = 60— - 432-^^ + 288-^^ . (51) 
" ^ 12r (12r)2 (12r)3 ^ ^ 

To first order in the density, the percolation threshold is then given by a straight line 

r = 2xiV ■ (52) 



^ = 7^ 10., , , on.. .2 • (^3) 



The next order already yields a more complex solution involving both X2 and Xs 

-6xiT^ + V Qr^^Wx3 + 6xfr - 1 8^2^ + 30xir^ 
X3 - 18X2^ + 30xir 
We then see that for r]_ < rj < ri_^_ with 

^ -Qjxl - 3x2) ± ySQixi - 3X2)^ - 120X1X3 .... 

60x? ' ^ ^ 

Eq. (l53l) has no real solutions. Clearly, the acceptable part of the solution is that for t] < t]^. 



B. CI integral equation 

The unphysical loss of the solution for the percolation threshold as obtained from 
the virial expansion is present also in the isotropic SHS model {6 = tt). This short- 
coming does not occur in an integral equation approach.— Within the Cn class of clo- 
sures a crucial role is played by the angular average of the cavity function at contact 
y = {y{ri2 = (T, r^i, ^2, ^i2)e^(l, 2))ni,t72,ni2- Its density expansion reads 

y^"" = yo + yiV + y2V^ + ■■■ , (55) 
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where Un = VniT) are related to the reduced virial coefficients 6„ by Eq. (HTl) . For i/q 
and yi they have already been computed in Eq. fl33l) . The percolation threshold is then 
given by 77^'^" = r where y*^" = (y'^"(ri2 = a, ili, ^2, ^^i2)e(l, 2))j7i,f72,!^i2 and ?/*^"(ri2 = 
cr, fii, f22, ^12) is the contact cavity function within the Cn approximation. 

Since = e then y = y and within the CO approximation [y'-'^ = yo) we find, 

T = xiV ^ (56) 
whereas within the CI approximation {y^^ = 1/0 + ViV) find, 



Now the loss of solution occurs between 



-3(X? - 4x2) ± V9(xf- 4X2)^ -120X1X3 
30x? 



^± = ^7^1:2 ■ (58) 



Note that at small values of 5 a gap may also appear in the CI percolation threshold for 
5 < 1.21 in the one-patch model and for 6 < 1.22 in the two-patch case. Fig. [T2] summarizes 
our findings and compares with MC simulations. From the figure we see that for 6 close to 
71 the percolation threshold of the two-patch model lies above that of the one-patch case at 
same total surface adhesive coverage, while the opposite trend is observed at lower 6. This 
mirrors our previous results for the coexistence curve. 

Another quantity which is useful to assess the onset of a phase transition is the average 
coordination number, defined by 

Z = p/^r,2(/.+ (l,2)>^^^^^_^^^. (59) 

One finds 

Z = 2^y^" , (60) 

which on the percolation threshold gives Z = 2. This prediction is compared with MC results 
in Fig. [12] where we show the average coordination number at the percolation threshold 
obtained from the MC simulations for the one- and two-patch models at 60% coverage. 

We are now in a position to summarize the phase diagram for one and two patches within 
the CI approximation. This is reported in Fig. [141 For two patches, the CI phase diagram 
(coexistence curve and percolation line) is compared with MC results in Fig. [15] both for the 
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full isotropic case (100% coverage) and for 60% coverage. Note that while the percolation 
line terminates at the point shown, the coexistence curve has a solution for the whole range 
of packing fraction considered. However, we have chosen to terminate the plot for the same 
value of 7] as the percolation line. 

VII. PHASE DIAGRAM AND ADDITION OF ADHESIVE BACKGROUND 

So far we have considered the case of adhesive patches on hard spheres. The disadvantage 
of this model is that there is no fiuid-fiuid transition below a certain surface coverage Xo- We 
have also shown that at fixed surface coverage the liquid more easily forms if the adhesion 
is distributed in different patches on the sphere surface for sufficiently large patches and we 
expect the opposite to be true for low coverage. 

One could argue that a more physical model should have a strongly directional potential 
mimicking e.g., active sites in a globular protein, in addition to an underlying isotropic 
attractive potential favoring a general fluid-fluid phase transition. To this aim, we modify 
our potential by adding a uniform adhesive background to each sphere on top of which a 
patchy potential of the type considered so far is active. This effect can be obtained by a 
simple substitution e — > 1 + Ae with A measuring the strength of adhesion on the patches, 
yielding 

Xi ^ 1 + Axi , (61) 
X2 ^ 1 + 2Xxi + X\2 , (62) 
X3 1 + 3Axi + 3X^X2 + A'X3 , (63) 

where 

X2 = (e(l,2)e(l,3))ni,Q2,n3lei2=vr/3 • (64) 

The phase diagram is now modifled as depicted in Fig. [16] where we have set A = 1 to be 
the strength of the patches throughout. In this case we see that even a small sticky patch 
(of amplitude 6 ~ 0.5) is sufficient to raise both the binodal and percolation threshold of the 
isotropic model. At equal coverage, the binodal and percolation threshold of the two-patch 
model lie below their one-patch counterparts, in agreement with the observed trend in the 
absence of background adhesion. 
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Note that the critical point is now less sensitive to the size of the patches because an 
isotropic SHS — rather than a hard sphere — is now the limiting 6^0. Indeed, the 

critical point does not move along t] while it covers the whole range < t < 2tc (A = 1), 
where Tc is the critical reduced temperature of the isotropic model (see Fig. [T7|) . The critical 
point shifts of the one- and two-patch models are now almost indistinguishable even though 
the crossing at 80% coverage still remains. 

VIII. CONCLUSIONS 

In this work we have studied, through integral equation theories and Monte Carlo sim- 
ulations, the structure, percolation and fluid-fluid coexistence curves of a model of hard 
spheres with one or two uniform sticky patches on their surface. Particles interact through 
an adhesive Baxter potential only if patches on different spheres are suitably oriented and 
as hard spheres otherwise. Unlike most previous studies, we have been able to analyze in 
some detail the dependence of the aforementioned quantities on the size of the patch and 
its interplay with the number of patches. 

The integral equation theory is based on the first two approximations (CO and CI) of a 
class of closures Cn which have already proved to provide a good qualitative representation 
of the exact behavior and are almost fully analytical. The comparison between the analytical 
work and Monte Carlo simulations indicates that CI yields a gratifying qualitative descrip- 
tion of the phase diagram notwithstanding the expected limitation due to its low density 
nature. While for the thermodynamics the results from the integral equation theories are 
exact within the given closure, for the percolation problem and the structure an additional 
orientational mean field approximation is necessary to decouple the orientational average. 

Radial distribution functions within the CO and CI integral equation theories exhibit a 
characteristic jump aX r = 2a (whose magnitude depends on the patch angle 5) and a cusp 
at r = 3(7. The coexistence and percolation lines move to lower temperature as the patch 
angle decreases from vr (the isotropic case) to (hard spheres). For a fixed surface coverage 
above approximately 80%, the curves of the two-patch case lie above the corresponding 
single-patch ones, while the opposite trend is observed below that point. We have suggested 
that this is due to two patches of sufficiently large size being able to form bonds to more 
particles than can a single patch. We have also argued that this reasoning does not apply 
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at low coverage, and that, in fact, the opposite situation might be expected. The crossover 
is not observed in the MC simulations within the range of adhesive coverage studied here 
(about 60%), but we cannot exclude the possibility for lower coverage, where the simulations 
converge very slowly. When an adhesive background is included in addition to the patches, 
both the liquid and percolating phase of the system are favored with respect to the isotropic 
case even in the presence of very small patches. 

In spite of the limited number of cases (one or two patches) addressed in the present work, 
our analysis suggests that both the total fraction of the surface covered by adhesion and the 
number of patches are crucial parameters in controlling the location of the critical point. In 
the limit of a single bond per patch, our analysis is consistent with a recent suggestion^ of 
a generalized law of corresponding states for non-isotropic patchy interactions. We remark 
that, from the purely theoretical point of view, there exist only few paradigmatic toy models 
with anisotropic interactions amenable to analytical or semi-analytical treatment. 

Our analysis can be regarded as complementary to recent investigations of the phase 
diagrams of globular proteins^'^'^ in that our starting point is the isotropic sticky hard 
sphere from which some adhesion is removed, rather that a hard sphere to which highly 
localized attractive spots are added. This approach goes beyond the limitation of one bond 
per patch, which is an essential feature of Wertheim thermodynamic perturbation theory. 
The price to pay is, of course, that only a qualitative agreement with MC simulations can 
be achieved. 

It would be interesting to extend the present work in some respects. In view of the diffi- 
culties of MC simulations in probing low coverage, a comparison with a numerical solution of 
a more robust closure such as, for instance, the Percus-Yevick approximation which has a full 
analytical description in the isotropic case, would provide a more quantitative assessment 
of the results presented here. Such a solution would also help to evaluate the (uncontrolled) 
angular decoupling approximation exploited in the present analysis of structure and the 
percolation threshold. Work along these lines is in progress and will be presented in a future 
publication. 
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APPENDIX A: THE LAW OF CORRESPONDING STATES 

Consider the simplest possible dependence A(si,S2) = Si ■ §2 and assume that e(l,2) = 
A(si,S2), i.e. the adhesion coefficient does not depend on ri2 = ri2/ri2. Within the Weeks- 
Chandler- Andersen perturbative expansion^ of the Helmholtz free energy A^^^ one finds 

/j(/lSHS_ A^S) ^ I ^(i)^(2)a«(ri,r2;r/)Ae(l,2) + 

J d{l)d{2)d{3)d{A) a(2)(ri, rs, rg, r/) Ae(l, 2) Ae(3, 4) + . . . (Al) 

where d{i) is a short-hand notation for dridQi, with dQi the average solid angle 
sin 6 idO i dip {/{Att), A^^ the Helmholtz free energy of the reference hard sphere (HS) sys- 
tem, 

Ae(l,2) = ^5(ri2-a) , (A2) 
i/r 

and the functions a*^") are expressed in terms of the correlation functions of the reference 
system which only depend on the packing fraction rj = irpa^/G, with p the density. We see 
then that the angular dependence in e factorizes and one finds 

i 

where Xi = / dQidQ2 2). This analysis shows how, in this case, the law of corresponding 
states holds. For example, if r = g{r]) is the spinodal or binodal of the SHS system with 
isotropic interaction (e = 1),— 1^1^ then the spinodal or binodal of the SHS with directional 
adhesion will be r = Xi9{v)y which will lie above that of the isotropic system if xi > 1 and 
below otherwise. 
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APPENDIX B: THE 3RD VIRIAL COEFFICIENT 

In this appendix we provide a derivation of Eqs. IQ, ([7]), ([8]), and ([9]). We start from the 
usual definition of the third virial coefficient 

B, = -^j dv, dr, dvs if (1, 2) / (1, 3) / (2, 3))^^ ..^^^^ (Bl) 

where, in fine with Eq. ([2]), the Mayer function can spht into two terms 

f{i,j) = fHsi^:J) + ^e{^,J)S{r,^-a). (B2) 

izr 

In the above equation we have set the HS part to the usual form 

/hs(^,j) = -e(a-r,,). (B3) 
Upon expanding the product, one can easily find 

AB3 = ^3 - ^3^^ = AB^^^ + ABf^ + ABf , (B4) 
where B^^ = Stt^ct^/IS is the HS result and 

ABi'^ = -i (^) J dr, dr2 dr, {Ms (1, 2) /hs (1, 3) e (2, 3) 5 (r^g - ^))n,,n,,ns (B5) 

A5f = -i(^)Y dr, dr, dr^ 

(/hs (1, 2) e (1, 3) 6 in, - a) e (2, 3) 6 {^3 - a)) (B6) 

(e (1, 2) 5 in, - a) e (1, 3) 5 (m - a) e (2, 3) 5 {m - ^))n,n„n3 • (BT) 

The above integrals are most conveniently evaluated in bipolar coordinates by introducing 
ri2 = r2 — ri and ri3 = — ri. This leads to r,^ = a/ + ^13 ^ 2ri2ri3ri ■ where it is 
most convenient to choose ri3 as the z axis. For AB^^ one finds 

AS$'^ = -{^) y"rfr23rfri2e(a-ri2)e(a-|ri2-r23|)5(r23-a)(e(2,3))^^_^|B8) 

Here one first performs the integration over ri2, which covers twice a spherical cap of height 
a/2 and then the straightforward integration over r23. Clearly the anisotropic part decouples, 
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thus yielding the isotropic part times Xi claimed. For ^ a little more care is necessary. 
One first obtains 

I dn^2dn^sS (^sJrl^ + rl,-2n,n,cose,2 - a) (e (1, 3) e (2, 3))^^^^^^^3 . (B9) 

After a first integration over ri3, an additional integration over cos 6*12 then requires ^13 < 7r/3 
(corresponding to the maximum available angle for all three particles in reciprocal contact). 
This also yields a normalization factor 4 = 1/ sin^(7r/6) in order to have the correct limit 
e{i,i) 1 for all («, j). The final result is 

Ai?f = / dfi,3 4 (e (1, 2) e (2, 3) ^ - ^13) ) , (BIO) 

thus yielding the isotropic part times X2 as reported in Eq. ([8]). An almost identical procedure 
also gives 

Ai?f = -(— ) ^-y^ dr,2rl,5 - a) j ^ rf(cosM 

5 (^7r?2 + ^2-2ri2acos^i2 - {e (1, 2) e (1, 3) 6 (2, 3))^^^^^^^^ (Bll) 

which, after an integration over the angular variables, leads to the desired decoupling for 
the anisotropic part Xs as given in Eq. (Q. Note that in this configuration all three spheres 
are necessarily touching and this fixes the angles % to a well defined value given in Eq. iQ. 
This completes the derivation of Eq. IQ. 



APPENDIX C: COEFFICIENTS X2 AND X3 FOR THE ONE-PATCH CASE 

Here we give the analytic expressions for the coefficients Qi and -Ri used in Eqs. f|T2l) and 
(fT3|) of the main text in terms of characteristic integrals which are then evaluated numerically. 
The basic procedure follows a similar analysis carried out in a different context,— which 
requires the calculation of the solid angle associated with the intersection of two identical 
patches on the same sphere as indicated in Fig. [181 For xs, one can easily see that for 6 < n/Q 
there is no possibility of intersection, even in the close-packed configuration. For 6 > 71/6 
the form of the resulting integral can be most conveniently written in slightly different ways 
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depending on the amplitude S of the patch. 



6 + 



lb. 

where the various terms are given in terms of the integrals 
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(i^ sin 6 arccos 
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sin 9 sin 27r /3 

cos(7r — 6) — cos 6* cos 7r/3 
sin ^ sin ir/3 



(CI) 

(C2) 
(C3) 
(C4) 



For example ab.i given in Eq. ( 1C2I) is the simplest integral resulting from the calculation of 
the overlapping region of the two cones of width 6 as depicted in Fig. [181 

For X2 an additional complication arises from the additional degree of freedom given by 
the fact that only two of the three spheres are (in general) in contact. One finds 
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APPENDIX D: COEFFICIENTS X2 AND xs FOR THE TWO-PATCH CASE 



Here we give the analytic expressions for the coefficients Q2 and R2 used in Eqs. (fTSll and 
IT6l) of the main text: 



(Dl) 
(D2) 
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LIST OF FIGURES 

Fig. 1 Summary of the vector notation used to define the model in the text. 

Fig. 2 Left paneh the adhesive patch model of Kern and Frenkel.— Right panel: patch surface 
coverage Xo ^ function of the patch angle 6. 

Fig. 3 Adhesion requires simultaneous alignment of patches (dark shading) on both spheres 
with the vector between their centers. The spheres on the left do not adhere, while 
those on the right do. 

Fig. 4 Configurations of three mutually bonded spheres, each possessing a large single patch 
(dark shading). The patch vectors point inward in the left panel and outward in the 
right panel. The latter case is only possible for 6 > Stt/G. Combinations of these 
arrangements are also possible. 

Fig. 5 Dependence of the coefficients Xi (^ = 1? 2, 3) on 6 for the one-patch model. 

Fig. 6 Configurations of three mutually bonded spheres, each possessing two patches (dark 
shading). On the left only one patch on each sphere is involved in the bonds; on the 
right both patches on each sphere are involved. The latter case is only possible for 
6 > Tf/3). Combinations of these arrangements are also possible. 

Fig. 7 Dependence of the coefficients Xi (^ = 1? 2, 3) on S for the two-patch model. 

Fig. 8 Radial distribution function for the one-patch model (n = 1) at 77 = 0.4 and r = 0.2 
within the CI approximation and for various values of the adhesive coverage nxo- 

Fig. 9 Comparison between the radial distribution function from Monte Carlo (MC) simula- 
tions and from the CO and CI approximations in the isotropic case {nxo = 100%, 
left-hand panel) and for an intermediate value of the single patch (n = 1) case 
with nxo = 80% (right-hand panel). Both sets of calculations were performed at 
r/(nxo)^ = 0.125 and pa^ = 0.35 corresponding to 77 = 0.183. . .. 

Fig. 10 Dependence on the adhesive coverage nxo of the binodal line calculated from the 
modified Carnahan-Starling free energy of Eq. (1371) (left-hand panel) and from the CI 
approximation of Eq. (l40l) (right-hand panel). The one-patch (n = 1) and two-patch 
(n = 2) systems are compared at the same total coverage. 
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Fig. 11 Dependence of the critical reduced temperature on the total adhesive coverage nxo for 
n = 1 and 2 patches, calculated from the CI approximation of Eq. PUI) and from MC. 
The inset shows the critical packing fraction. 

Fig. 12 Percolation thresholds for various total adhesive coverages as calculated from the CI 
approximation (lines), Eq. (!57|) . and MC simulation (points). The one-patch {n = 1) 
and two-patch (n = 2) cases are compared at the same total coverage. 

Fig. 13 Average coordination numbers at the percolation threshold for the one- and two-patch 
models at 60% coverage, obtained through MC. The continuous line is the prediction 
from the integral equation theory. The isotropic case is also reported for comparison. 

Fig. 14 Phase diagram in the CI approximation, for various values of the adhesive coverage 
nxo- The one-patch (n = 1) and two-patch (n = 2) models are compared at the same 
total coverage. 

Fig. 15 Comparison of the CI approximation with MC simulation (dots) for the phase diagram 



of particles with two patches. The MC isotropic phase diagram is taken from Ref. 
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Fig. 16 Dependence of the binodal line of patchy adhesive spheres with a background adhesion 
on the total surface coverage nxo of the patches: on the left as calculated from the 
modified Carnahan-Starling free energy of Eq. (|37l) : on the right as calculated from 
the CI approximation, Eq. fHUjl . 

Fig. 17 Dependence of the critical reduced temperature on the adhesive coverage in the patchy 
model with a uniform adhesive background, calculated from the CI approximation, 
Eq. fl4UI) . The inset shows the behavior of the critical packing fraction. 

Fig. 18 Basic geometry for the calculation of Ri{6) in Eq. fICll) . The required solid angle is 
the overlap of the two cones of width 6 (the darkly shaded region in the sketch). 
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